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Abstract 

We present a lattice gas cellular automaton (LGCA) to study spatial and temporal 
dynamics of an epidemic of SIR (susceptible-infected-removed) type. The automaton 
is fully discrete, i.e. space, time and number of individuals are discrete variables. The 
automaton can be applied to study spread of epidemics in both human and animal 
populations. We investigate effects of spatial inhomogeneities in initial distribution 
of infected and vaccinated populations on the dynamics of epidemic of SIR type. We 
discuss vaccination strategies which differ only in spatial distribution of vaccinated 
individuals. Also, we derive an approximate, mean-field type description of the au- 
tomaton, and discuss differences between the mean-field dynamics and the results of 
LGCA simulation. 

1. Introduction 

Since the publication of Kermack and McKendrick epidemic model |l|] in 1927, mathemat- 
ical epidemiology developed an extensive body of literature. A vast majority of proposed 
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models, to study dynamics of epidemics, is based on ordinary differential equations. These 
models assume that the population mixing is strong, hence concentrations of effected types 
of population (for example, susceptible, infected, removed) are spatially homogeneous. Thus, 
these models neglect spatial aspects of the epidemic process. Models which rely on partial 
differential equations, such as [0], abandon the assumption of homogeneous mixing and al- 
low to study the geographical spread of epidemics, yet they still pose some serious problems. 
They treat the population as continuous entity, and neglect the fact that populations are 
composed of single interacting individuals. This can lead to very unrealistic results, such as, 
for example, endemic patterns relaying on very small densities of individuals, named in P| 
"atto- foxes" and "nano-hawks" . 

Models based on interacting particle systems eliminate these shortcomings of traditional 
methodologies. They treat individuals in biological populations as discrete entities, allow to 
introduce local stochasticity, and are well suited for computer simulations. Several models 
of this type have been suggested in a recent decade, including stochastic interacting particle 
models 0], and models based on cellular automata or coupled map lattices p, |^, |^, ||, |^. 

In this article, we present lattice gas cellular automaton (LGCA) for SIR (susceptible- 
infected-removed) type of an epidemic. Our aim is to demonstrate that the LGCA approach 
provides an interesting and potentially fruitful alternative to other methods. In real biologi- 
cal populations, both animal and human, contact processes among infectious and susceptible 
individuals and their movement play a vital role in spreads of epidemics of an infectious dis- 
ease. Infectious diseases spread because infectious and susceptible individuals mix together. 
They move, meet each other and through a contact process they transmit an infection. 
Hence, among other factors the spread of infectious diseases strongly depends on patterns of 
mobility in populations. In LGCA, population mixing arises directly from motion of individ- 
uals and susceptible individuals can become infected only if they meet infectious individuals. 
Additionally, a LGCA allows to investigate effects of spatial inhomogeneities in population 
concentrations on the dynamics of epidemic processes and vaccination strategies. We will 
demonstrate examples of such effects. In this article we do not aim to construct a LGCA 
of a particular disease but rather we want to study generic features of LGCA methodology 
and its suitability in the context of epidemiology. To illustrate the ideas and discuss fur- 
ther possible developments we selected SIR (susceptible-infected-removed) epidemic type for 
which we constructed a LGCA. We derived approximate mean-field type description of the 
automaton and simulated and analyzed the automaton dynamics. We compared the pre- 
dictions obtained from the mean-field description with those obtained from the automaton 
simulations. Additionally, we used automaton to study various vaccination strategies. 



2. Individual-based SIR model on a lattice 

We will construct a LGCA for an epidemic of SIR type, for which we can assume that 
a population consists of three types of moving and interacting individuals, of type r G 
{S,I,R}, susceptible, infected, and recovered. The proposed automaton is a special case a 
lattice gas cellular automaton for reaction-diffusion systems, described in detail in [TT 



We tile the physical space, in which an epidemic takes place, by regular hexagonal cells 
h{r) with centers at discrete space variables r. We chose the hexagonal cells in order to 
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Figure 1: The array of hexagonal cells with a dashed lattice superimposed on it. The lattice 
depicts the channels through which the individuals are transported from cell to cell. The 
figure is reproduced from [|12 . 



avoid spurious invariants in the dynamics of the automaton, also known as parity problem 



10| , pj]] . However, other types of tiling, such as square cells, can be used as well. We assume 



that the distance between centers of adjacent cells is 1. Let 

Ci =< cos(i — 1)11/3, sin(i — l)7r/3 > (1) 

be a unit vector, for each i = 1, 6. If we connect the center of every hexagon h{r) with the 
centers of the neighbouring hexagons r + Cj, where i = 1, 6, then we obtain a hexagonal 
lattice structure Ch, with the lattice coordination number m = 6. In the case of square cells 
the lattice coordination number would be m = 4. Figure |I] (reproduced from [Q) shows that 
the centers of the hexagons h{r) become the nodes of the hexagonal lattice Ch- The cells 
represent some area of the physical space in which the individuals mixed, interact and can 
move from one place (cell) to another. The individuals can move among the cells only along 
the channels corresponding to edges of the lattice Ch- Except where confusion might arise, 
we will identify a cell h{r) with its center r, hence with the node r of the lattice Ch- Hence, 
we can say that the individuals reside at the nodes of the lattice and move along its edges 
(channels). While extensions to arbitrary boundary conditions are straightforward, we limit 
our considerations to periodic boundary conditions. 

At the initial time, k = 0, particles representing individuals 5", / and R are distributed 
randomly and independently, among the cells h{r), r G Ch, according to probabilities given 
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by their initial concentrations. The initial distribution of particles is such that there are 
at most six particles, regardless of the individual type, per cell h{r) (four when the cells 
are squares). The LGCA discrete dynamics is constructed in such a way that the total 
number of individuals at a given node is restricted to stay between and m (m = 6 or 4, 
respectively, for hexagonal or square lattice) per each cell during the time evolution of the 
automaton. Since a single lattice node represents a small region of space, this means that 
m can be interpreted as the carrying capacity of that region. In other words, the number of 
individuals that this region can support cannot exceed m. 

The time evolution of the automaton takes place at discrete time steps. At each time step 
k, an evolution operator £ is applied, simultaneously and independently of the past, to all 
lattice nodes. The evolution operator S governs the dynamics of the automaton, which arises 
from sequential applications of the three basic operations contact C, randomization TZ and 
propagation V. Hence, the evolution operator can be written in terms of these operations as 
the superposition 

E = VonoC. (2) 

Each of the operations C, TZ and V captures some aspects of the epidemic process and their 
actions are described as follows. 

1. As a result of an application of the contact operation C, individuals can change their 
type, meaning that susceptible individuals can become infected, and infected individu- 
als can recover. More precisely, each susceptible individual at a node r, independently 
of other individuals, can become infected with probability 1 — (1 — r)^' , where uj is a 
number of infected individuals at the node r, and r G [0, 1]. Similarly, each infected in- 
dividual at the node r, independently of other individuals, can recover with probability 
a, where a G [0, 1]. 

2. As a result of using the randomization operation TZ, applied at each node independently 
of the other nodes, a population of individuals residing at a node r is randomly re- 
distributed among edges/channels originating from the node r. Through the selected 
channels, allowing at most one individual per channel, individuals will move in the 
propagation step from the node r to the neighbouring nodes. The process of redis- 
tribution of individuals is purely probabilistic one and it contributes to modeling the 
mixing process of individuals. 

3. In the propagation step, governed by the operator V, individuals simultaneously move 
from their nodes to the neighbouring ones through the channels assigned to them in 
the randomization step. The movement of individuals is purely deterministic in the 
propagation step. 

The rationale behind the probability of becoming infected in the first step being 1 — (1 — 
r)"^ can be explained as follows. We assume that, at a given node, a susceptible individual 
contacts all infected individuals at that node and that all infected individuals are infectious 
one. If the probability of infection per contact is r, then the probability of not getting 
infected after contacting each of uj infected individuals is (1 — r)"^. The probability of 
getting infected, therefore, is 1 — (1 — r)"-^. 
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Note that the mechanism of contracting an infection described here imphes that the in- 
cubation period is short enough to be neghgible, meaning that a susceptible who contracts 
the disease at time k becomes infectious immediately, and can infect others at time k + 1. 
While this is a convenient assumption, it will have to be relaxed in more realistic models, or 
in models tailored for a specific disease. Also, note that the duration of the disease (number 
of time steps from contracting an infection to recovery) has geometrical distribution with a 
mean value of 1/a time steps. Again, in realistic disease-specific models a different distri- 
bution will have to be adopted, depending on the particular disease. As we said earlier the 
spread of infectious diseases strongly depends on mixing patterns in the populations. In the 
LGCA that is being presented here, mixing is of diffusive type. It arises from randomiza- 
tion of directions of motion of individuals in the randomization step TZ and the movement 
of individuals in the propagation step V. In the more realistic models the movement of 
individuals will have to be appropriately modified and additionally, we will have to include 
inflow, outflow, birth and death processes. 



3. Mean-field approximation 

In order to gain some insight into the dynamics of the LGCA defined in the previous sec- 
tion, we will now proceed to construct approximate equations describing the automaton 



dynamics. Using the formalism and methodology introduced in [^], it is possible to write 
compact microdynamical equations corresponding to the evolution operator S. Moreover, 
by making appropriate approximations, it is possible to derive the LGCA discrete lattice- 
Boltzmann equations and the corresponding partial differential equations describing the 
automaton macroscopic dynamics. Since rigorous derivation of these equations is rather in- 



volved, we refer the interested reader to |Tl[]. We take here a less rigorous, but more intuitive 
approach. 

Let us assume that the total number of nodes in the lattice Ch is A^, and at time k there 
are Ns{k) susceptible, Nj{k) infected, and Nfi{k) recovered individuals on the entire lattice. 
Since we do not allow for inflow, outflow, birth and death processes, the total population, 
denoted by Ntot, remains constant in time. Define ^^-(r, /c) to be a number of individuals 
of type r at a node r at time k. If we assume that the automaton dynamics is spatially 
homogeneous (well stirred), then the expected value of nr{r, k) is independent of r and equals 

{nrir,k)) = ^. (3) 

Under the same assumptions, regardless at which cell the susceptible individual is located, 
the expected number of infected individuals at his cell is Nj/N, thus the probability that he 
becomes infected is 

1 _ (1 _ r fi(k)/N_ 

Hence, the expected number of susceptible individuals who become infected in a single time 
step is equal to 

A^5(A;)[1 - (1 -r)^^(^)/^]. (5) 
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Similarly, the expected number of individuals who become recovered in a single time step 
is aNj{k). When the population is well stirred this yields that the expected number of 
individuals of each type r G {S, I, R} at time k + 1 is 



Ns{k + 1) 
Nj{k + 1) 

NR{k+l) 



iV5(A:)(l-r)^^('^)/^, 

Nj{k) + Ns{k)[l - (1 - rf'O^y^] - aNj{k), 
Nn{k)+aNi{k). 



(6) 



For small r, taking Taylor expansion 



+ ..., 



(7) 



N 2m 



and by keeping only the first two terms and defining Prik) = 



Nr{k)/N, we obtain 



Ps{k + 1) 
pi{k + l) 
pRik + 1) 



ps{k) - rps{k)pi{k), 
pi{k) + rps{k)pi{k) - 
PR{k) + api{k). 



api{k), 
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The above equations are quite similar in their structure to the ordinary differential equations 
obtainable from the classic Kermack-McKendrick modelQ, as described, for example, in [p^j . 
The similarity lies in the fact that the gain in the class of infected individuals occurs at a 
rate proportional to the density of infectives and susceptibles, in analogy to the mass action 
law in chemical kinetics. Of course, this is valid only for small values of r and under the 
assumption of strong mixing, which is not always realistic. 

In reality, the epidemic process spreads quite differently than the mean-field approxima- 
tion predicts. Figure ^ compares the number of infected individuals as a function of time 
as observed in the LGCA simulations with the mean- field approximation (|^). The simula- 
tions have been performed on a hexagonal lattice with lO'^ nodes, using r = 0.3, a = 0.2. 
The initial configuration consisted of 16000 susceptibles and 100 infected individuals for each 
simulation. At the initial time for each simulation the individuals were randomly and uni- 
formly distributed on the lattice. The simulation curve in Figure |^ represents average over 
50 experiments. The mean-field approximation (MFA) predicts much faster initial spread 
than observed in the LGCA simulations. This can be easily explained by realizing that 
the MFA assumes perfect mixing, meaning that an infected individuals can always infect 
some susceptibles residing at their node. In simulations, since the mixing is limited, it often 
happens that the number of infected individuals in the vicinity of an infected individuals is 
larger than average and, similarly, the number of susceptibles is smaller than average. Or, 
in other words, infected individuals are often grouped together in small regions of space, 
while in other spatial regions there are no infected individuals at all. Therefore, the effective 
force of infection is smaller than predicted by MFA. We should also note that the fraction of 
susceptibles who eventually became infected is 67% for simulations and 94% for mean-field, 
meaning the epidemic is less severe when the mixing is weaker, even though it lasts longer, 
as can be seen in Figure @. 

^The Kermack-McKendrick model is based on integral equation with general time-kernel for infectivity. 
System of ODE we are referring to can be obtained from this integral equation as a special case. 
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Figure 2: Number of infected individuals as a function of time, comparison of LGCA sim- 
ulation results with the mean-field approximation. The simulation curve represents average 
of 50 LGCA experiments. Each LGCA simulation has been performed on a hexagonal lat- 
tice with 10^ nodes with r = 0.3, a = 0.2, and with 16000 susceptibles and 100 infected 
individuals uniformly and randomly distributed on the lattice at = 0. 



4. Effects of spatial distribution 

In order to illustrate the importance of spatial distribution of individuals in spread of epi- 
demics, let us consider the following problem. Let us assume that in some region of space, to 
be denoted by A, several cases of an infectious disease have been reported. In order to limit 
the spread of the disease, we would like to vaccinate the population. However, we have only 
M doses of a vaccine at our disposal, where M is less than the number of susceptible indi- 
viduals in the population. The natural question arises how should these doses be distributed 
in the population to minimize the severity of the epidemic? To simplify the problem, we 
assume that the vaccine is immediately effective, therefore vaccinated individuals become 
members of the recovered group immediately after vaccination. Additionally, we assume 
that the region A is a circle of radius 20 on a hexagonal lattice of 10^ nodes. 

We compare two alternative vaccination strategies. In the first one, to be referred to 
as a "uniform strategy", we vaccinate M individuals selected from the entire population 
of susceptibles randomly and independently of each other. In the second strategy, to be 
called a "barrier strategy", we vaccinate all individuals in a ring surrounding the region 
A. The thickness of the ring is selected in such a way that the total number of vaccinated 
individuals to be M. In all subsequent discussions we assume M = 1000. Figures |]a and 
^ (page 13) show initial distribution of individuals in both these strategies, "uniform" and 
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Figure 3: Number of infected individuals as a function of time, under "uniform" and "barrier" 
vaccination strategies. Each simulation curve represents average of 100 LGCA experiments. 
For each type of vaccination strategy each simulation has been performed on a hexagonal 
lattice with 10^ nodes with r = 0.3, a = 0.2, and with 16000 susceptibles uniformly and 
randomly distributed on the lattice aX k = 0, and 10 infectives uniformly and randomly 
distributed in the circle of radius 20 on the lattice at /c = 0. 

"barrier", respectively. Note that the number of individuals in S, I, and R classes is the same 
in both cases, so that the spatial distribution of individuals is the only difference between 
configuration in Figures ^ and The color coding used is green for S, red for J, and 
black for R. Consecutive snapshots of the LGCA dynamics in Fig. ^ reveal that the time 
evolution of these two initial configurations is quite different. One immediately notices that 
the final configuration in the case of the "uniform" vaccination strategy, after 1000 time steps 
(Fig. ^), contains many more blue pixels than the corresponding final configuration of the 
"barrier" vaccination strategy (Fig. ^i). This means that many more susceptible individuals 
contracted the disease in the course of the epidemic when the "uniform" vaccination strategy 
was used than in the case of the "barrier" vaccination strategy. To be more precise, out of 
initial 16000 susceptible individuals out of whom 1000 have been vaccinated at the start of 
epidemic, 52% became infected when the vaccination was "uniform", and only 12% when 
the "barrier" vaccination strategy was employed. This can be seen from Figure ^ where 
the number of infected individuals is plotted as a function of time step, averaged over 100 
LGCA experiments. 
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Figure 4: Average severity of epidemics of an infectious disease {Nji{oo) — Nji{0)/Ntot), 
calculated as average over severities of outbreaks of epidemics of the disease, as a function of 
barrier permeability /. For each / the average (■) has been taken over 100 LGCA simulations 
of outbreaks of epidemics, each starting with identical initial configuration, and r = 0.3 and 
a = 0.2. 

5. Barrier's permeability and severity of the epidemics 

Since the barrier appears to be a very effective strategy, one could ask how permeability of the 
barrier influences the dynamics of the epidemics. In what follows, we will vary permeability 
of the barrier by vaccinating not all susceptibles in the ring surrounding the region A, but 
only a fraction / of randomly selected susceptibles in the ring. Thus, if / = 0, there is no 
single vaccinated individual in the ring, while / = 1 corresponds to the "barrier" described 
in the previous section, where all susceptibles in the ring are vaccinated. 

We define the severity of an outbreak of an epidemic of an infectious disease as (A''r(oo) — 
Nii{0))/Ntot, where Nii{oo) denotes a total number of removed when a total number of 
infected in LGCA simulation reaches a steady state and Nfi{0) denotes the total number of 
removed (vaccinated individuals) at the start of an epidemic (in an initial configuration in our 
case). Since in our model we do not allow for birth, death and migration of individuals then 
the number of removed will remain constant from the moment when the number of infected 
reaches zero. Let Ns{oo) be a total number of susceptibles at the end of an outbreak of the 
epidemic. Since Nm = Ns{0) + Ni{0) + Nn{0) = Ns{oo) + Nr{oo), then Nr{oo) - Nr{0) = 
-^s(O) + Nj{0) — Ns{oo) gives us the total number of susceptibles who got infected in the 
course of the epidemic. Hence, the average severity of an epidemic of infectious disease can 
be measured by taking the average ((A''r(oo) — NR{0))/Ntot) of the severities of outbreaks of 
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specific epidemics of tlie infectious disease. It tells us on average what chance a susceptible 
has to become infected in a course of an epidemic. The severity of an epidemic or of an 
outbreak of the specific epidemic of an infectious disease depends on many factors including 
the infectivity probability per contact r, the recovery probability a and spatial distribution 
of infected and removed at the outbreak of an epidemic. 

In order to examine how the average severity of an epidemic of an infectious disease 
depends on the barrier permeability we performed a series of simulations varying /. For 
each / we measured the average severity of an epidemic as the average over severities of 
100 epidemic outbreaks each with the infectivity probability per contact r = 0.3 and the 
recovery probability a = 0.2 and strating from the initial conditions described as follows. 
For each barrier permeability / each simulation has been performed on a hexagonal lattice 
with 10^ nodes, and with 10 infectives uniformly and randomly distributed in the circle A 
of radius 20 on the lattice at = 0, and with 16000 susceptibles uniformly and randomly 
distributed on the lattice, from which in the ring surrounding the circle and consisting of 
1000 susceptibles we "vaccinated", at = 0, randomly and independently a fraction / of 
them. Results of this LGCA experiment, which show the average severity of an epidemic 
of SIR type as a function of barrier permeability / are presented in Figure ^ It is rather 
remarkable that as / approaches zero, the average severity of an epidemic of SIR type comes 
close to 49%, not much different than the corresponding value for the "uniform" vaccination 
strategy, when r = 0.3 and a = 0.2. This means that the "uniform" vaccination strategy is 
practically not better than no vaccination at all in this case. Even quite sparse "barrier" is 
much more effective than uniform random vaccination. 

6. Conclusions 

We presented a lattice gas cellular automaton for studying spreads of epidemics of SIR type. 
We derived an approximate mean-field type description of the automaton, and discussed 
differences between the mean-field approximation and the results of the simulation using 
LGCA. We also investigated what effects can have spatial inhomogeneities in the distribution 
of various types of populations on the dynamics of the epidemic process. We demonstrated 
that the severity of an epidemic can strongly depend of an initial spatial distribution of 
vaccinated individuals. 

The lattice gas cellular automaton discussed in this work is certainly simplistic. It does 
not take into account many complexities which have to be considered when one attempts to 
construct a more realistic automaton model for an epidemic. However, the important feature 
of our model is explicitness of mixing and contact processes. Unlike models based on partial 
differential equations, our model is individual-based, and the spread of the infection occurs 
due to the motion of individuals and their interactions. In such models it is quite straight- 
forward to introduce different, non-diffusive type of motion, and investigate the effects of 
resulting mixing on the dynamics of the epidemic process. For example, we are currently 
investigating periodic motion with some amount of randomness, that might better represent 
the behavior of individuals in human populations. Results of such experiments will be pre- 
sented elsewhere. Here we would like only to reiterate that the description of non-diffusive 
motion with partial differential equations is usually very difficult or impossible, excluding 
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trivial cases such as linear transport. Individual-based models are much more suitable for 
this purpose, and we hope that they eventually will help to built epidemic spatial models, 
tailored for specific diseases, with high degree of realism. 
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Figure 5: Snapshots of LGCA dynamics evolution under "uniform" (left) and "barrier" 
(right) vaccination strategies. 
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